Fast Color Space Transformations Using Minimax Approximations 
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Abstract 

Color space transformations are frequently used in image processing, graphics, and visualization applications. In many 
cases, these transformations are complex nonlinear functions, which prohibits their use in time-critical applications. In this 
paper, we present a new approach called Minimax Approximations for Color-space Transformations (MACT). We demonstrate 
MACT on three commonly used color space transformations. Extensive experiments on a large and diverse image set and 
comparisons with well-known multidimensional lookup table interpolation methods show that MACT achieves an excellent 
balance among four criteria: ease of implementation, memory usage, accuracy, and computational speed. 



Introduction 
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OSolor space transformations are commonly used in various image processing, graphics, and visualization tasks for decoupling 
^5aminance and chromaticity information, ensuring approximate perceptual uniformity, or achieving invariance to different 
. ijiiaging conditions such as viewing direction, illumination intensity, and highlights. Since color devices are usually provided 
^^ith direct RGB signal input and output, the RGB color space is generally the source in these transformations. 
k>< _ In many cases, color space transformations are complex nonlinear functions, which prohibits their use in time-critical 
^Applications. In this paper, we present a new approach called Minimax Approximations for Color-space Transformations 
C^MACT). We demonstrate MACT on transformations between the RGB space and three popular color spaces: CIELAB, HSI, 
and SCT. There are various reasons behind the choice of these particular spaces. First, all three spaces decouple luminance 
and chromaticity information, which makes them suitable for tasks including enhancement [1] [5] , rendering [3] , noise removal 
in E], segmentation [61 [71 [8], and object recognition [9l [10]. Second, an important feature of CIELAB is its approximate 
perceptual uniformity, an essential requirement for expressing color differences in applications such as color quantization 
mesh denoising [Tl|, and maximum contrast color set design [T3]. Third, HSI and SCT are designed to match human intuition, 
which makes them useful for manual color selection |14] . Last, the first two components of HSI are shown to be invariant to 
viewing direction, surface orientation, illumination direction, and illumination intensity |10] . 

The rest of the paper is organized as follows. Section [5] gives the related work. Section [3] presents the use of minimax 
approximation theory in speeding up color space transformations and the experimental results. Finally, Section [4| gives the 
conclusions. 



2 Related Work 

Traditionally, color space transformations in digital imaging systems have been implemented using lookup tables (LUTs) that 
require some form of multidimensional interpolation. The most commonly used 3D LUT interpolation methods are trilinear, 
prism, pyramidal, and tetrahedral (15j . These methods will be explained in ^3.5.31 In this section, wc briefly review the 
interpolation methods that appear less frequently in the literature as well as some alternative approaches. 

Chang et al. [T^] developed a method called Sequential Linear Interpolation (SLI) that uses a partially separable grid 
structure, which allows the allocation of more grid points to the regions where the function to be interpolated is more nonlinear. 
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This scheme often resuhs in more accurate transformations at the expense of increased computational cost when compared to 
trihnear interpolation. Kacker et al. |17j proposed a wavelet based method that uses a multiscale grid structure. This method 
is shown to achieve lower maximum error but higher average error when compared to SLI. Gupta and Gray |18j presented a 
method called Maximum Entropy Estimation, which is a generalization of tetrahedral interpolation. This method is shown to 
be more accurate than tetrahedral interpolation in general while being at best only half as fast. Hemingway |19| described a 
method based on n-simplexes that is slightly faster than tetrahedral interpolation. However, the accuracy of the method is not 
discussed in the paper. 

Neural networks have been applied to the color space transformation problem in a number of studies [501 [HI HH US] • They 
have the advantages of being more flexible and requiring less memory. However, they require training and parameter tuning, 
and are prone to overfitting [TB]. 



3 Approximate Color Space Transformations 
3.1 Overview of Minimax Approximation Theory 

Given a function /, we would like to approximate it by another function g so that the error (e) between them over a given 
interval is arbitrarily small. The existence of such approximations is stated by the following theorem: 

Theorem 3.1 (Weicrstrass) Let f be a continuous real-valued function defined on [a,b], i.e. f G C[a,6]. Then Ve > there 
exists a polynomial P such that j|/ — P|| < s, i.e. Vx £ [a, 6], — Pix)\ < e. 

This is commonly known as the minimax approximation to a function. It differs from other methods, e.g. least squares 
approximations, in that it minimizes the maximum error (e) rather than the average error: 

e = max \f{x) - P{x)\ (1) 

xG[a^b] 

A similar theorem establishes the existence of a rational variant of this method |24| . Let n > be a natural number and let 

Pn{[a,b]) = {ao + aix + . . . + Gnx'^ : X e [a,b], Ui E R, i = 0,l,...,n} (2) 

be the set of all polynomials of degree less than or equal to n. The set of irreducible rational functions, {[a, b]), is defined 
as: 

^;^(K&]) = {^:P(2;) q(x)eF„([a,6])| (3) 

where p and q have no common factors. Then |24] : 

Theorem 3.2 For each function f e C[a, b], there exists at least one best rational approximation from the class R^^-^ ([a, 6]). 

This theorem states the existence of a rational approximation r* € ([a, b\) to a function / G C[a, b] that is optimal in 
the Chebyshev sense: 

max |/(a:)-r*(x)| =dist(/,i?:;,) (4) 

where dist (/, i?"J denotes the distance between / and ([a, b]) with respect to some norm, in our case, the Chebyshev 
(maximum) norm. Regarding the choice between a polynomial and a rational approximant, it can be said that certain functions 
can be approximated more accurately by rationals than by polynomials. Jean-Michel Muller explains this phenomenon as 
follows "It seems quite difficult to predict if a given function will be much better approximated by rational functions than by 
polynomials. It makes sense to think that functions that have a behavior that is 'highly nonpolynomial' (finite limits at ±oo, 
poles, infinite derivatives, . . .) will be poorly approximated by polynomials." |25| . 

In this study, the Remcz Exchange Algorithm is used to calculate the minimax approximations. The reader is referred to 
|24[ 125] for more information on the theory of minimax approximations and [2 6) for the implementation details of the Remez 
algorithm. 



3.2 CIELAB Color Space 
3.2.1 Color Space Description 

CIELAB is an approximately uniform color space standardized by the CIE (Commission Internationale de I'Eclairage) in 1976 
[27] . The L* component represents the lightness, whereas a* and b* represent the chromaticity. The transformation between 
the RGB and CIELAB color spaces is comprised of three steps: (1) conversion from nonlinear R'G'B' to linear RGB, (2) 
conversion from linear RGB to CIEXYZ, and (3) conversion from CIEXYZ to CIELAB. 
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In order to convert the nonlinear R'G'B' values to the linear RGB ones, inverse gamma correction (ITU-R BT.709) is 
performed: 



^ \ 255' 255' 255 J 
k e {R, G, B} = 




< fc' < 0.081 
0.081 < fc' < 1 

The conversion from linear RGB to CIEXYZ (ITU-R BT.709) is given by: 

X = 0.412391i? + 0.357584G + 0.180481S 
Y = 0.212639i? + 0.715169G + 0.072192B 
Z = 0.019331i? + 0.119195G + 0.950532B 

Finally, the conversion from CIEXYZ to CIELAB is given by: 

L* = 116,f{Y/Yo) - 16 

a* = 500 [/ (X/Xo) - f (Y/Yo)] 

b* = 200[f{Y/Yo)-f{Z/Zo)] 

' t^l^ t > 0.008856 

^ 7.787t + 16/116 t < 0.008856 

Here, Xg, Yq^ and Zq arc the tristimulus values of the reference white. For the illuminant D65 these are: 



(5) 



(6) 



(7) 



fit) 



Zo 



: 0.950456 
1.0 

1.089058 



The distance between two pixels x and y in the CIELAB space is given by: 

DciELAB{^,y) = [Ll - L*Y + (a* - a*)^ 
3.2.2 Approximation of the Cube-Root Function 



(8) 



(9) 



^ and (jni) can be implemented efficiently using LUTs as commonly seen in the literature [TS]. The cube- root function (cbrt) 
in ([7]) is the main factor that influences the computational time of the transformation. The probabilities of calling this function 
can be calculated from an image that contains every possible color in the 24-bit RGB space (see ij3.5.ip : 



P{X/Xo > 0.008856) = 0.999626 
P iY/Yo > 0.008856) = 0.999165 
p \z/Zo > 0.008856) = 0.997184 



(10) 



These high values suggest that we can accelerate the transformation substantially if we can devise a fast approximation for 
the cube-root function (see Figure [T]). Table [T] shows the coefficients of the minimax polynomials of various degrees. Here, n 
and e^aa: represent the degree of the polynomial and error of the minimax approximation, respectively. 



0.0 0.2 0.4 0.6 0.8 1.0 



Figure 1: Cube-root function in the interval [0.008856, 1] 
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Tabic 1: Minimax polynomials for the cube-root function 



n 




aa 


ai 


a2 


as 


04 


a 5 


2 


1.271154C-01 


1.268979C-01 


2.393873 


-1.647669 








3 


9.787829C-02 


9.787826C-02 


4.057495 


-7.388864 


4.331370 






4 


8. 1111500-02 


8. 1111330-02 


5.926004 


-2. 0171650+01 


2. 8330700+01 


-1. 3247280+01 




5 


7. 0029560-02 


7. 0029100-02 


7.961214 


-4. 3293520+01 


1. 0631820+02 


-1. 1356850+02 


4. 3582680+01 



Tabic 2: Minimax rationals for the cube-root function 



n 


m 






ai 


02 


^3 


a4 


2 


2 


2. 0609960-03 


6. 3096550-03 
4.482646O-02 


5. 7857820-01 
1.175862 


1.591005 
9.596879O-01 






2 


3 


7. 2102310-04 


2. 5007050-03 
1. 9787010-02 


3. 4471130-01 
8. 5427970-01 


1.942708 
1.664540 


-2. 5032670-01 




3 


2 


5.931593e-04 


1. 7765190-03 
1. 4322560-02 


2. 6323230-01 
6. 7799980-01 


1.751297 
1.706260 


3. 8367090-01 




2 


4 


3.107735e-04 


1. 3178990-03 
1.124254O-02 


2. 3901130-01 
6. 7266790-01 


2.099395 
2.184656 


-7. 5111500-01 


2.229717e-01 


3 


3 


1.858694e-04 


4. 3708890-04 
3.912364O-03 


9. 5269520-02 
2.954084O-01 


1.252009 
1.717143 


1.302733 
6. 3434080-01 




4 


2 


2. 3346880-04 


7. 5893020-04 
6. 6447230-03 


1. 5197840-01 
4.506424O-01 


1.663584 
2.030622 


8. 3680750-01 


-1. 6572690-01 


3 


4 


8. 5398630-05 


1. 6836670-04 
1. 6108640-03 


4. 6676750-02 
1. 6179740-01 


9. 1068120-01 
1.494070 


1.810577 
1.218468 


-1. 0794510-01 


4 


3 


8. 0529200-05 


1. 3496730-04 
1. 2992500-03 


3. 8320790-02 
1. 3454200-01 


7. 8701740-01 
1.330358 


1.799062 
1.365369 


2. 0711700-01 


4 


4 


3. 8569300-05 


3. 9272830-05 
4.0221000-04 


1. 3923180-02 
5. 4145360-02 


4. 1147390-01 
8. 2215260-01 


1.734853 
1.800167 


8. 6792230-01 
3.513617O-01 



It can be seen that the error values are quite high and as the approximation degree is increased, the accuracy doesn't 
improve significantly. This suggests that rational functions might be better suited for this approximation task. Table [2] shows 
the coefficients of the minimax rationals of various degrees. Here, each pair of adjacent rows corresponds to a rational function 
of a particular degree (n, m) in which the first and second rows represent the numerator and denominator, respectively. It can 
be seen that minimax rationals can accurately represent the cube-root function. 



3.3 HSI Color Space 
3.3.1 Color Space Description 

HSI (Hue-Saturation-Intensity) is an intuitive alternative to the RGB space |27| . It uses approximately cylindrical coordinates, 
and is a non-linear deformation of the RGB color cube. The hue H is a function of the angle in the polar coordinate system 
and describes a pure color. The saturation S is proportional to radial distance and denotes the purity of a color. Finally, the 
intensity I is the distance along the axis perpendicular to the polar coordinate plane and represents the brightness. 
The transformation between RGB and HSI is given by: 



H 



0.5{R-G+R-B) 



(11) 



y^{R-Gf + {R-B){G-B) 

a {B>G) ' 

5 = 1-3 mm{R,G,B)/{R + G + B) 
I = {R + G + B)/3 

where arccos denotes the inverse cosine function. Kender PH] proposed a fast version of ([TT|) that gives numerically identical 
results. This transformation involves fewer multiplications and no square root operation: 



if i? > S and G> B 

H = ^ + arctan 
else if G > i? 

H = n + arctan 
else if B > G 

arctan 
else \i R> B 

else 

H = undefined 



V3(G-R) 
G-B+R-B 



V3{B~G) 



B-R+G~R 



V3[R-B) 



R~G+B-G 



(12) 
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Tabic 3: Minimax polynomials for the inverse tangent function 



n 




oo 


ai 


a2 


as 


04 


^5 


1 


6.907910C-03 


5.959793C-03 


1.782975 


7.497879C-01 








3 


3.654156C-03 


-3.654076C-03 


1.884080 


-9. 8055830-01 


1. 4305800-01 






4 


1.286371C-03 


-1.286369C-03 


1.796716 


-4. 9589690-01 


-6. 9274040-01 


4. 4215410-01 




5 


1.801311C-04 


-1.801283C-04 


1.739333 


-2. 0398480-02 


-2.065512 


2.052837 


-6.591729C-01 



where arctan denotes the inverse tangent function. The distance between two pixels x and y in the HSI space is given by: 



Dhsi i^,y) = ys'i + Sy - 2sxSy cos6' + (ix - iy)'^ 

\h^-hy\ if \h^-hy\ < IT (13) 

277 — \h^ — hy I otherwise 



3.3.2 Approximation of the Inverse Tangent Function 

We decided to approximate Kender's transformation (fT2)) rather than the original one pT|) due to several reasons. First, these 
two give identical results. Second, (fT2|) is computationally cheaper than (fTTj) . Third, as will be seen in the next subsection, 
the inverse tangent function in (|12|) is easier to approximate when compared to the inverse cosine function in (jlip . 

In (fT2|) . all of the cases involve multiplication by the constant v^- Therefore, a multiplication operation can be avoided 
by approximating arctan (a/S^;) ■ Note that the argument of this function can also be negative. This can be handled using the 
following identity: 

arctan (x) =— arctan (— a;) (14) 

In (|12p . the inverse tangent function receives its arguments from the interval [—1.0,1.00 Figure [2] shows a plot of the 
function in the second half of this interval. 



0.6- 



0.2- 



0.0- 




0.0 0.2 4 0.6 0.8 1.0 

X 



Figure 2: Inverse tangent function in the interval [0, 1] 

Owing to its highly linear behavior, this function can be accurately approximated by low-order polynomials. Table [3] shows 
the coefficients of the minimax polynomials of various degrees (n). 

3.4 Spherical Coordinate Transform (SCT) Color Space 
3.4.1 Color Space Description 

The Spherical Coordinate Transform (SCT) is defined as |29j : 

L = VR^ + G^ + B^ 

ZA = arccos (B/L) ^]^5j 
ZB = arccos ( ^,4^^) ) 

where L represents the luminance, and angles ZA and ZB represent the chromaticity. 
The expression for ZB can be simplified using trigonometric manipulations: 

ZB = arctan (G/R) (16) 



^This holds prior to the multiphcation with \/3. 
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This formulation is computationally advantageous in that it avoids a multiplication and a sine operation. 

Unfortunately, it is not easy to define a perceptual distance function in SCT. Although various formulae have been developed 
to calculate the distance between two points lying on the same spherical surface, these cannot be used in SCT. This is because 
pixels of different brightness in this space lie on different spherical shells. Therefore, we decided to use the following alternative 
approach. Given two pixels in SCT, the inverse transformation (jl7|) is applied to switch back to the RGB space. The pixels are 
then converted to CIELAB so that the distance between them can be calculated using ([9]). This indirect method of distance 
calculation is likely to introduce additional errors. However, minimax approximations for the elementary functions involved in 
(jlSp and (jl6p can be devised to obtain an arbitrarily accurate approximate transformation. 

i? = L sin (ZA) cos (ZB) 

G = Lsin(ZA)sin(ZB) (17) 
B = L cos (ZA) 

3.4.2 Approximation of the Inverse Cosine Function 

In the expression for ZA, the inverse cosine function receives its arguments from the interval [0, 1]. Figure [3] shows a plot of 
this function. 



1.6- 




0.0 0.2 4 0.6 0.8 1,0 

X 



Figure 3: Inverse cosine function in the interval [0, 1] 

Unfortunately, approximating the inverse cosine function in this interval is not easy because of its behavior near 1. This 
can be circumvented using the following numerically more stable identity for x > 0.5: 



arccos(a;) = 2 arcsin I a/ 0.5(1 — x) ) (18) 




Figure 4: Inverse sine function in the interval [0,0.5] 

In (I18|) . the inverse sine function (arcsin) receives its arguments from the interval [0,0.5]. Figure 3] shows a plot of this 
function. In order to avoid two multiplication operations, the following function can be approximated instead: 



arccos(.T) = 2 arcsin (y/-\/2) 
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Tabic 4: Minimax polynomials for the inverse sine function 



n 




Of) 


ai 


a2 


as 


04 


(15 


4 


2.097814C-05 


2. 0977970-05 


1.412840 


1. 4298810-02 


6. 7043610-02 


6.909677C-02 




5 


2. 3705400-06 


-2.370048C-06 


1.414434 


-3.300037C-03 


1. 3546700-01 


-3. 9942590-02 


6. 0995020-02 


Table 5; Minimax polynomials for the inverse cosine function 


rn 


A 


ao 


ai 


a 2 


as 


a4 


as 


4 


1.048949C-05 


1.570786 


-9.990285C-01 


-1. 4298990-02 


-9. 4813350-02 


-1.381942C-01 




5 


1.186403C-06 


1.570798 


-1.000156 


3.299810C-03 


-1. 9157800-01 


7. 9882310-02 


-1.725177C-01 



In (fT9|) . the argument y falls into the interval [O, l/V2]- Table |4] shows the coefficients of the minimax polynomials of 
various degrees (n). 

On the other hand, it can be seen from Figure [3] that the inverse cosine function is highly linear in the interval [0, 0.5] and 
can be accurately approximated by polynomials. Table [5] shows the coefficients of the minimax polynomials of various degrees 
(m). 

3.4.3 Approximation of the Inverse Tangent Function 

In (jl6p . the inverse tangent function receives its arguments from the interval [0,255]. Figure [5] shows a plot of this function. 
Note that the degenerate cases, i.e. R = G = (ZB = undefined) and R = 0, G > {ZB = 7r/2), are not reflected in the plot. 



Figure 5: Inverse tangent function in the interval [0, 255] 

It can be seen that the function exhibits high variability in this large interval. In order to obtain an accurate low-order 
polynomial approximation, the domain can be divided into two as follows: 



X = G/R 
arctan(a;) 



f{x) = arctan(.T) ii G < R 

gix) = 7r/2 — arctan (l/a:) otherwise 



(20) 



The second part of follows from: 



arctan(a;) = 7r/2 — arctan(l/a;) for a; > 



(21) 



When G < R, the function / receives its arguments from the interval [0,254/255]. On the other hand, when R < G the 
inverse tangent function within g receives its arguments from the interval [1/255,1]. Tables [6] and [7] show the coefficients of 
the minimax polynomials of various degrees (r) for the functions / and g, respectively. 

As in the case of ([TOl) . the probabilities of calling the inverse sine and inverse cosine functions can be calculated from an 
image that contains every possible color in the 24-bit RGB space: 



-Parosin = P {B / L > 0.5) = 0.557723 
P — 1 — P — 442277 

arccos ^ arcsm u.^'iz.z. i i 

In contrast, the inverse tangent function in (jl6p is almost always called (except for when R 
ZB = undefined): 



Pa, 



= 1-P{R = G = 0) = 1 - 



256 

2563 



= 0.999985 



(22) 

G = 0, in which case 
(23) 



Note that ZB when i? = and G > 0. 
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Tabic 6: Minimax polynomials for the function / 



r 




ao 


ai 


a2 


£13 


04 


tJ5 


4 


1.036515C-04 


-1. 0365080-04 


1.003740 


-1.773538C-02 


-3.390563C-01 


1.386796C-01 




5 


2.073939C-05 


2.073866G-05 


9.982666C-01 


2.352573C-02 


-4.506862G-01 


2.635050G-01 


-4. 9208220-02 



Table 7: Minimax polynomials for the function g 



r 




ao 


ai 


02 


aa 


04 


05 


4 


1. 0516430-04 


1.570917 


-1.004004 


1. 8856940-02 


3. 3731590-01 


-1. 3779300-01 




5 


2. 0121040-05 


1.570769 


-9.981253C-01 


-2.440212C-02 


4. 5289210-01 


-2. 6591810-01 


5. 0162280-02 



3.5 Experimental Results 

3.5.1 Calculation of the Color Space Transformation Errors 

In order to calculate the accuracy of the presented approximate color space transformations, we used a 4, 096 x 4, 096 image 
(henceforth referred to as RGBWMillion]^, which contains 16, 777, 216 unique colors, i.e. every possible color in the 24-bit 
RGB space. Tables [8IIT0I show the average execution time (in seconds jl over 1,000 identical runs and the average (s'^yg) and 
maximum (Ema^.) transformation errors for CIELAB ([9]), HSI (|13p . and SCT ^ & pT]) . respectively. Note that the rows of 
Table [8] are sorted on [n + m), since the execution time is proportional to the total degree of the polynomials in the rational. 
On the other hand, the rows of Table [TOl are sorted on the degree of the inverse tangent approximation (r). This is because 
the probability of calling this function ([23| is much higher than that of the other two elementary functions ([22]). Therefore, 
the execution time is mainly influenced by the degree of the inverse tangent approximation. 

We also performed the exact transformations on RGB16Million and calculated the average execution time over 100 runs. 
The results were 60.331s, 13.485s, and 27.599s for CIELAB, HSI, and SCT, respectively. Comparing these values with those 
given in Tables [8lll01 we can see that the proposed approximations provide substantial computational savings. 

3.5.2 Calculation of the Computational Gain Values 

It can be seen from Tables [Sl IlOl that only marginal computational gains can be obtained using lower-order approximations 
without significantly compromising the accuracy of the transformation. Therefore, in this subsection, we consider only the 
highest-order approximations for each transformation. However, lower-order approximations might be preferable depending on 
the application requirements. 

In order to calculate the computational gain for each approximate transformation, a set of 100 high quality RGB images 
was collected from the Internet. The set includes images of people, animals, plants, buildings, aerial maps, man-made objects, 
natural scenery, paintings, sketches, as well as scientific, biomedical, synthetic images and test images commonly used in the 
literature. 

On each image, the exact transformations were performed 100 times and the average execution times were calculated. The 
same was done for the corresponding approximate transformations with 1,000 runs. The computational gain is calculated as 
the ratio of the average execution times. Table [TT] shows the statistics over the entire image set. It can be seen that the 
computational gain observed in RGB16Million also applies to the case of real- world images. 

3.5.3 Comparison with 3D Lookup Table Interpolation Methods 

In this subsection we compare MACT with the most commonly used LUT interpolation methods, namely trilinear, prism, 
pyramidal, and tetrahedral. These methods involve three steps: packing, extraction, and interpolation |15j . Packing is a 

^Available at [http: //brucelindbloom. com/downioads/RGB16Million. tif .zip 

•^Programming language: C, Compiler: gcc 3.4.4, CPU: Intel Pentium D 2.66Ghz 



Table 8: Comparison of CIELAB approximations 



n 


TO 


Avg. Time 


^ avq 




2 


2 


1.031872 


0.787741 


2.225399 


2 


3 


1.134568 


0.295577 


0.774404 


3 


2 


1.114776 


0.239033 


0.636929 


2 


4 


1.221390 


0.131411 


0.331015 


3 


3 


1.165900 


0.061127 


0.193565 


4 


2 


1.197828 


0.090846 


0.247291 


3 


4 


1.291526 


0.018092 


0.085046 


4 


3 


1.291570 


0.014927 


0.079935 


4 


4 


1.292872 


0.003201 


0.036481 
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Table 9: Comparison of HSI approximations 



n 


Avg. iimc 


s 

avq 


s 

max 


2 


0.865328 


0.002523 


0.007134 


3 


0.865575 


0.001259 


0.003693 


4 


0.866254 


0.000448 


0.001320 


5 


0.873251 


0.000063 


0.000190 



Table 10: Comparison of SCT approximations 



n 


m 


r 


Avg. Time 


^avq 


^max 


4 


4 


4 


1.330927 


0.006596 


0.024273 


4 


5 


4 


1.371811 


0.006553 


0.024273 


5 


4 


4 


1.377924 


0.006266 


0.022820 


5 


5 


4 


1.419315 


0.006224 


0.021810 


4 


4 


5 


1.427337 


0.002139 


0.008371 


4 


5 


5 


1.462627 


0.002011 


0.008371 


5 


4 


5 


1.471486 


0.001381 


0.005793 


5 


5 


5 


1.507812 


0.001254 


0.004543 



process that divides the domain of the source space and populates it with sample points to build the LUT. The extraction step 
aims at finding the location of the input pixel and extracting the color values of the nearest lattice points. The last step is 3D 
interpolation, in which the input point and the extracted lattice points are used to calculate the destination color specifications. 
Essentially, 3D interpolation is a repeated application of linear interpolation [15) : 

a; = (1 — a)xo + axi (24) 

where Xq and Xi are the spatial coordinates of the two known points and a is the interpolation coefficient. The 3D inter- 
polation is the step in which the aforementioned LUT interpolation methods differ. Trilinear interpolation uses 8 neighboring 
lattices, whereas prism, pyramidal, and tetrahedral interpolations use 6, 5, and 4 neighbors, respectively. The greater the num- 
ber of neighboring lattices used in a method, the higher the computational requirements and accuracy. Since the formulations 
of these methods are mathematically involved, the interested reader is referred to the relevant literature jl5j . 

We implemented two versions of each LUT method: standard and caching. The former is a direct implementation of the 
mathematical formulation, whereas the latter pre-computcs the differences between neighboring lattices and stores these in 
each node. Table [T2] shows the average execution times on RGB16Million over 1,000 runs and the transformation errors for 
CIELAB. It can be seen that even though these methods range from being 1.32 to 2.87 times faster than MACT, their accuracy 
is 18 to 161 times lower. In addition, these methods require extra storage of up to 10.7 MBs. Note that LUT size only affects 
the accuracy of the transformation and the storage requirements and, in theory, it should not affect the computational time. 
However, it can be observed from the table that as the LUT size is increased, the computational time increases as well. This 
is most likely due to the limited size of the processor cache. 

As for HSI and SCT, since they include angular components (H component in HSI, /LA and /LB components in SCT), 3D 
interpolation in these color spaces involves the interpolation of angular (circular) data given by |30j : 

e = arctan ^ ^ - ^) ^in (^o) + « sin (^,) \ ^^^^ 



(1 — a) cos (6*0) + a cos (6*1) 

where Oq and 9i arc the two known angles and a is the interpolation coefficient. Due to the trigonometric and inverse- 
trigonometric functions in (|25|) . it turns out that 3D interpolation in these spaces is computationally more expensive than the 
original transformations given in (|lip and (|T5|) . For example, trilinear interpolation in HSI takes about 84s on RGB16Million 
(eLg = 0.000777, = 1-086255). 

Note that even if the angular data was quantized into, say 360 steps, prccomputing (|25p would not be possible since the 
value of the interpolation coefficient a is not known a priori. In addition, the inverse tangent function cannot be approximated 



Table 11: Computational gain statistics for each transformation 



Transformation 


Min 


Con 

Max 


iputational ( 
Mean 


3ain 
Stdev 


Median 


CIELAB 


8.761736 


40.629135 


35.575373 


4.608590 


37.192729 


HSI 


4.511637 


16.460248 


11.437106 


2.439482 


11.669715 


SCT 


4.493543 


22.232897 


15.782453 


2.110820 


16.162434 
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Table 12: Comparison of 3D LUT interpolation methods 



Method 


LUT 
Size 


Avg. Time 
(standard) 


Avg. Time 
(caching) 


7 

e 

avg 


7 ' 

e 

max 


Trilinear 


93 
173 
333 


0.901673 
0.908753 
0.982880 


0.619681 
0.630465 
0.716958 


0.359677 
0.100026 
0.025089 


5.409771 
1.987040 
0.653506 


Prism 


93 
173 
333 


0.856652 
0.860880 
0.931349 


0.609452 
0.614681 
0.677319 


0.349128 
0.095150 
0.024172 


5.409771 
2.343905 
0.826651 


Pyramidal 


93 
173 
333 


0.850409 
0.872813 
0.940143 


0.540481 
0.566271 
0.637049 


0.319163 
0.088116 
0.022466 


5.872591 
2.909481 
1.042665 


Tetrahedral 


93 
173 
333 


0.689341 
0.704121 
0.773102 


0.450413 
0.486654 
0.561076 


0.284339 
0.077870 
0.020253 


5.783921 
2.702448 
1.185788 



using a minimax polynomial since its argument is not bounded. In short, the expensive interpolation operation offsets the 
computational advantage of 3D LUT interpolation in angular color spaces. 

4 Conclusions 

In this paper we proposed MACT, a novel approach to speed up color space transformations based on minimax approximations. 
Advantages of MACT include ease of implementation, negligible memory requirements, extremely good accuracy, and very 
high computational speed. Comparisons with commonly used 3D LUT interpolation methods revealed that MACT yields 
significantly more accurate transformations at the expense of slightly higher computational requirements. Although MACT 
was applied to three particular color space transformations, it can easily be adapted to other transformations, e.g. RGB to 
CIELUV transformation, that involve computationally expensive mathematical functions. 

Implementations of the fast color space transformations described in this paper will be made publicly available at |http : / /j 
[sourcef orge .net/projects/f ourier-ipal. 
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